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1-10 keV spectrum of AXPs/SGRs forms as the thermal photons emitted 
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^ . The anomafous X-ray pulsars and soft 7-repeaters are peculiar high- 

energy sources believed to host a magnetar, an ultra-magnetized neutron star 

I with surface magnetic field in the PetaGauss range. Their persistent, soft X- 

, ray emission exhibit a two component spectrum, usually modeled by the su- 

\^ . perposition of a blackbody and a power-law tail. It has been suggested that the 

(N 

' by the cooling star surface traverse the magnetosphere. Magnetar magneto- 

00 ; 

, spheres are, in fact, likely different from those of ordinary radio-pulsars, since 

the external magnetic field may acquire a toroidal component as a consequence 
• 1—1 . 

' of the deformation of the star crust induced by the super-strong interior field. 

H ' 

' In a twisted magnetosphere, the supporting currents can provide a large opti- 

cal depth to resonant cyclotron scattering. The thermal spectrum emitted by 
the star surface will be then distorted because primary photons gain energy 
in the repeated scatterings with the flowing charges, and this may provide 
a natural explanation for the observed spectra. In this paper we present 3D 
Monte Carlo simulations of photon propagation in a twisted magnetosphere. 
Our model is based on a simplified treatment of the charge carriers veloc- 
ity distribution which however accounts for the particle collective motion, in 
addition to the thermal one. Present treatment is restricted to conservative 
(Thomson) scattering in the electron rest frame. The code, nonetheless, is 
completely general and inclusion of the relativistic QED resonant cross sec- 
tion, which is required in the modeling of the hard 20-200 keV) spectral 
tails observed in the magnetar candidates, is under way. The properties of 
emerging spectra have been assessed under different conditions, by exploring 
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the model parameter space, including effects arising from the viewing geom- 
etry. Monte Carlo runs have been collected into a spectral archive which has 
been then implemented in the X-ray fitting package XSPEC. Two tabulated 
XSPEC spectral models, with and without viewing angles, have been pro- 
duced and applied to the O.f-f keV XMM-Newton EPlC-pn spectrum of the 
AXP CXOU J f 647-4552. 

Key words: Radiation mechanisms: non-thermal - stars: neutron - X-rays: 
stars. 



1 INTRODUCTION 



Over the last few years, increasing observational evidence has gathered in favour of the ex- 
istence of "magnet ars" , i.e. neutron stars (NSs) endowed with an ultra-strong magnetic 
field {B 10^^ — 10^^ G), much higher than the critical threshold at which quantum 
electro-dynamical (QED) effects become important {Bc-it ~ 4.4 x 10^^ G). The ex i stence 
of t hese objects has been f irst p roposed in the early '90s by iDuncan &: Thompson! (119921 ) 



and 



Thompson fc DuncanI (119931 ). who suggested that, soon after the core collapse follow- 



ing the supernova explosion, convective motions can strongly amplify the seed magnetic 
field via helical dynamo action. The magnetar model, initially developed to describe the 
phenomenology of the so-called soft 7-ray repeaters (SGRs), namely the emission of strong 
bursts, the fast spin period evolution and the persistent X-ray luminosity, is currently be- 
lieved to successfully reproduce the properties of another class of pecuhar NSs, the anoma- 
lous X-ray pulsars (AXPs). Although alternative models, invoki ng accretion from a fossil 
disc, are not completely ruled out by obse rvations as yet (see e.g. IVan Paradijs et al.lll995 : 



Chattereiee et al.l 



2000 



five AXPs (Gavriil et al. 



Krimm et al 



*erna et al. 



2OO2I : 



20001) . the recent detection o 



" SGR-like bursts from 



Kaspi et al. 



2003 



Woods et al. 



2005 



Kaspi fc Gavriil 



2006 



20061 ) has strengthened the connection between the two groups and pushed 



forward the interpretation of AXPs as magnetars. 

Both classes of sources, SGRs and AXPs, are characterized by spin periods in a narrow 
range (5-12 s), a typical persistent X-ray luminosity of ~ 10^^-10^^ ergs~^, no evidence 
for Doppler shifts in the light curve, lack of bright optical companions (favouring an inter- 



pretation in terms of isolated objects), and a spin-down in the range 10 ^^-10 ss ^. In 
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particular, the magnetar scenario appears promising in providing an alternative mechanism 
(namely the ultra-strong magnetic field) to power their high X-ray luminosity, which can 
not be otherwise explained in terms of more conventional processes, as accretion from a 
binary companion or injection of rotational energy in the pulsar wind/magnetosphere. Be- 
sides, measurements of period and period derivative, assuming that spin-down is associated 
to magneto-dipolar losses, are strongly suggestive of the presence of an ultra-strong magnetic 
field, B > B„a- 

The soft X-ray spectra of AXPs are generally well described by a two component model, 
consisting of a blackbody with kT ~ 0.4-0.5 keV, and a power-law with photon index 2 < 



r < 4 fe. 



Woods &: Thompson! I2OO6I . and references therein). In some cases, SGRs spectra 



have been fit with a single power-law component, but recent deep o bservations showed that , 



also for these sources, a blackbody component is often required ( iMereghetti et al. 



2005b . 



20061 ). Despite the fact that the blackbody plus power-law model has been routinely applied 



to magnetar candidate spectra for several years, attempts to provide a physical interpretation 
for these two components have just begun. 

Recently it has been proposed that this phenomenological spectral model mimics a sit- 
uation in which soft seed photons emerging, for instance, from the neutron star surface are 
boosted to higher energies by efficient resonant cyclotron scattering (RCS) from magneto- 
spheric charged particles, l eading to the formation of a power-law hi gh-energy tail. The basic 
idea has been discussed by Thompson. Lyutikov &: Kulkarni (j2002 ) (TLK in the following), 
who suggested that a possible difference between SGRs/AXPs and standard radio-pulsars 
is that in the former the internal magnetic field is highly twisted, up to ~ 10 times the 
external dipole. Stresses imparted to the star crust by the strong toroidal component of the 
internal magnetic field cause the crust to deform. This produces, in turn, a displacement 
of the footpoints of the external magnetic field lines with the net result that, at intervals, 
the external (initially dipolar) field may acquire a toroidal component, i.e. it may twist 
up as well. Twisted magnetospheres are threaded by currents, substantially in excess of the 
Goldreich- Julian current. As shown by TLK, charge carriers may provide large optical depth 
to resonant cyclotron scattering so that soft (thermal) photons produced at the star surface 
gain energy through repeated collisions with the moving charges. Since the electron distri- 
bution is spatially extended, and the resonant cross section depends on the local value of B, 
it is expected that repeated scatterings lead to the formation of a high-energy tail, instead 
of a narrow cyclotron line. At least qualitatively, this scenario may also explain the correla- 
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tion between spectral hardening, luminosity and increase in bursting/ glitc hing activity that 



has been recently discoyered in the long-term evolution of a few sources (iMereghetti et al. 



2nn5b 



Rea et al 



2005 



Campana et al.l 120071 ). 



The problem of computing the X-ray spectrum emerging from twiste d magnetospheres 



has be en previously tackled using a simplified one- dimensional approach by 
(120061 ). a nd a systematic a pplication to X-ray data has been presented by 



vutikov fc Gavriil 



Rea et al 



(see also 
sented by 



. lea et al 



mm 

2007a[) . More r ecent ly, 3-D Monte Carlo calculations have been pre- 



Fernandez fc Thompson! (120071 ). although these spectra have never been applied 



to fit X-ray observations. Both these investigations treat resonant cyclotron scattering in 
the non-relativistic regime, and neglect electron recoil, i.e. use the resonant cross section in 
the particle frame in the (magnetic) Thomson limit. 



Interestingly, thanks to INTEGRAL, it has been recently found that AX^ 
exhibit very hard high-energy tails (F ^ 1) which can extend up to ^ 200 keV 




den Hartog 



Gotz et al. 



et al. 



2006 



Revnivtsev et al 



2004 



Mereghetti et al 



2005al : 



's and SGRs 



Kuiper et al 



Molkov et al 



20061 ). This discovery come somewhat as a surprise, being the persistent 



spectra of these sources below ~ 10 keV rather soft, and changed our view of magnetars, 
suggesting that their luminosity might well be dominated by the hard, rather than soft, X- 
ray component. The origin of such high energy tails is presently unclear, but, again, most of 
the scenarios proposed so f a r invo ke emission from magnetospheric particles. Quite recently. 



Thompson fc Belobodorovl (120051 ) discussed how soft gamma-rays may be produced in a 



twisted magnetosphere, suggesting two different mechanisms: either thermal bremsstrahlung 
emission from the surface region heated by returning currents, or synchrotron emission from 
pairs created higher up (~ 100 km) in the magnetosphere. However, an alternative possibility 
is that the high-energy tails are again created by resonant magnetic Compton up-scattering 
of soft X-ray photons. In order to boost efficiently soft photons up to a few hundred keVs, 



scattering must occur off a non-thermal 



bly located close to the stellar surface (IBaring fc Harding 



population of relativist ic electrons (or pairs), possi- 



20071 ). A quantitative calculation 



of the expected spectra, which necessarily requires a correct description of relativistic effects, 
has not been put forward as yet. 

In this paper, the first in a series devoted to investigate the X-/soft 7-ray persistent 
spectrum of magnetar candidates, we lay out the physical bases of our model and present a 
Monte Carlo code which is used to follow the spectral modifications as the soft seed photons 
get progressively up-scattered in the magnetosphere of an ultra-magnetized neutron star. Our 
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present goal is to test, by direct comparison with observations, if RCS spectra are capable of 
accounting for the observed properties of the soft X-ray emission (< 10 keV) of SGRs/AXPs. 
To this end, we adopt a non-relativistic (Thomson) description for the scattering process. 
However, the numerical scheme is completely general and is explicitly designed to incorporate 
the fully QED cross sections and to deal with more complex magnetic configurations. The 
former, together with an application to the hard X-ray tails detected by INTEGRAL, will be 
the scope of forthcoming papers (Nobili, Turolla & Zane in prepar ation). In many respects 
the p resent investigation follows an approach similar to that of [Fernandez Sz Thompson 
(120071 ). and we will refer to this paper for some useful expressions. The two treatments, 
however, differ in a number of ways. In particular, the present model includes the angular 
and frequency dependence of seed photons in a more general way and a different prescription 
for the current velocity distributions. Differences and similarities between the two methods 
will be discussed along the paper, when relevant. 

The paper is organized as follows. In ^ we lay out and scrutinize the physical bases of 
our model. The Monte Carlo method and its coding is described in ^ while in §l]we present 
the computed spectra and discuss their properties. The implementation in XSPEC of our 
model is described in [|5] where also a preliminary fit is reported. Discussion follows in ^ 



2 THE MODEL 

In this section we discuss in some detail the main ingredients used in our computation of 
the soft (~ 0.1-10 keV) X-ray spectrum emitted by magnetar candidates. 

2.1 External magnetic field geometry 

The first ingredient of our computation is a prescription for the magnetic field geometry. 
Monte Carlo techniques are suitable for handling complicated 3D configurations, and our 
code is completely general from this point of view. However, for the sake of simplicity, in this 
paper we restrict ourselves to the axially symmetric twisted magnetosphere configurations 
studied by TLK, in which case the (numerical) solution of the magnetostatic, force-free 
equilibrium is straightforward. Accordingly, we report here only those expressions that are 
needed to facilitate the reading of this paper, and refer to TLK for all details. 

The starting point is the force-free equation j x B = where j and B are the current 
and the external field respectively. Under the assumption of axial symmetry, this equation 
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can be written as V x i? = a{V)B with V = V{r, 9) the flux parameter. A major simphfi- 
cation arises by restricting to self-similar configurations, V = Vor~'^F{cos6), in which case 
the problem reduces to the solution of a second order eigenvalue differential equation for 
F{cos6), that can be solved numerically for each value of the parameter ^ p ^ 1. The 
latter univocally fixes the magnetic configuration, a part for a scale factor Bpoie (see below) . 
The boundary conditions are chosen in such a way that the resulting axially-symmetric con- 
figuration corresponds to a core-centered, twisted, dipolar field (see ^for a discussion). The 
polar components of the magnetic field are then (see again TLK for all details) 

^ Bpoie f Rns\^^^ dF 



Bo 



Bpole 


f Rns 


2 


\ r 


Bpole f 


Rns\ 


2 V 


r J 



B^ = B, 



C 



pF_ 
sin^ 



p{l+p)_ 

where the constant C is an eigenvalue which depends on p only, Rjsis is the neutron star 
radius and -Bpoie is the value of the magnetic field at the pole. The net twist angle is defined 

as 

A^„_.= lim2r|*^ (2) 

00^0 Jeo Bg sm6 

and is a function of the parameter p. As a consequence, either p or A(^n-s can be used to 
label each model in the sequence. 

2.2 Magnetospheric currents 

Once the magnetic structure is known, in the force-free approximation the spatial density 
of the magnetospheric particles is automatically fixed by 

where (/?) is the average charge velocity (in units of c; see below). The above expression 
gives the co-rotation charge density of the space charge-limited flow of ions and electrons 
from the NS surface, that, due to the presence of closed loops in a twisted field, is much 
larger than the Goldreich- Julian density, ricj- Moreover, it is important to note that, while 
a space charge-limited flow with n = ncj requires currents flowing in opposite directions 
from the two poles, for the case at hand there is a well defined flow direction which is the 
same from north to south. This breaks the symmetry between the two star hemispheres, 
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and implies that the observed spectrum will be different when viewed from the north or the 
south pole. Clearly, because of charge neutrality, the electron current must be balanced by 
ions flowing in the opposite direction. However, ions are heavier, they are not lifted much 
in the magnetosphere and tend to move closer to the star surface. Photons may scatter 
off ions, but this is likely to giv e rise at most to a narrow ab sorption feature at the ion 
cyclotron energy (see TLK and [Fernandez Sz Thompson! 120071 ). For this reason, the ion 
current is not considered here, together with pair creation, that can further complicate the 
relation between charge and current density by introducing bidirectional flows (see ^ for 
a discussion). In a genuinely static twist {dB/dt = 0) the electric and magnetic fields are 
orthogonal. This implies that the voltage drop between the footpoints of a field line vanishes 
since =0, so that there is no force that can extract particles from the surface and lift 
them against gravity thus initi ating the current ?'r = cV x B/4 n requested to support the 
twist. However, as discussed in iBeloborodov &: ThompsonI (120071 ) . once implanted, the twist 
has necessary to decay precisely to provide the potential drop required to accelerate charges. 
A non-vanishing E\\ is maintained by self-induction and the twist evolution is regulated by 
the balance between the conduction current j and Jb, dE\\/dt = 47r(jB ~ j)- If J < Jb the 
magnetosphere becomes charge starved and -Ey grows at the expenses of the magnetic field, 
injecting more charges into the magnetosphere. On the other hand, when j > Jb the field 
decreases reducing the current. The magnetosphere is then in dynamical (quasi)equilibrium 
with ?' ^ ?'r over a timescale < tdprm,, where tdecay ~ a few years is the twist decay time 
(IBeloborodov fc ThompsonI 120071 ). 

The second key ingredient is the velocity distribution of the magnetospheric charges. This 
is a crucial and still largely unexplored issue (see however IBeloborodov fc ThompsonI 120071 . 
and Nevertheless, in a strong magnetic field the electron distribution is expected to be 
largely anisotropic: e~ stream freely along the field lines, while they are confined in a set of 
cylindrical Landau levels in the plane perpendicular to B. In order to mimic such scenario, we 
assume a 1-D Maxwellian distribution at a given temperature Tg, superimposed to a bulk 
motion with velocity Vbuik, as measured in the stellar frame. The (invariant) distribution 
function turns out to be 

drie He exp{-'y'/Qe 



^e/e(r,7/5) 



(4) 



d{^(3) 2iri(i/ee) 

where 7' = •yybuiki'^ — PPbuik), ©e = kTe/nieC^, Ki is the modified Bessel Function of the 
first order and fe = 'y~^n~^dne/d(3 is the momentum distribution function. We consider 
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Te and [3i,uik = v^uik/c as free parameters in our model, and, although this is definitely a 
simplification, we assume that both do not de pend on position. This expression differs from 
that used by [Fernandez fc Thompson! (120071 ) (their equation [19]) inasmuch they do not 
include the effects of collective (bulk) velocity (which is necessary to reproduce the current 
flow), but only those of the e~ local velocity distribution (either thermal, as in the present 
case, or non-thermal). In other words, we assume that electrons move isothermically along 
the field lines but, at the same time, they receive the same boost from the electric field. 
Even in the lack of any detailed information about the charge accelerating mechanisms, we 
consider our choice more realistic. 



2.3 Scattering cross sections 



Scattering off free electrons in the presence of a strong magnetic field has been extensively 
treated in the literature. The non-relativistic {B ^ Bcrit) expressions for the scatteririg cross 
sections in the Thomson limit (i.e. neglecting electron recoil) were derived by I Ventura! (119791 ) 



(see also 



Meszaros 



19921). The comple t e QED co mpton cross sections have been presented by 



Heroldl (119791 ). iDaugherty fc Hardingl (Il986l ). and lHarding &: Daughertyl (Il99ll ). The scatter- 



ing cross section depends on the incident photon polarization state and, in general, it must 
be computed by summing over the (infinite) virtual intermediate Landau states. Moreover, 
proper account has to be made for the electron spin transition and for the possibility that 
scattering leaves the electron in an arbit rary excited state (Raman s cattering). This leads 
to quite cumbersome expressions (see e.g lHarding &: Daughertylll99ll ). even if one restricts 
to the resonant part of the completely differential cross section. On the other hand, under 
the typical conditions expected in a twisted magnetosphere, soft photons {huj ~ 1 keV) 
will undergo resonant scattering when uj ^ ujb and this happens only where the field has 
decayed to a value B ~ 10^^ G ^ Bcrit- Electron recoil starts to be important when the 
photon energy in the electron rest frame becomes comparable to the electron rest energy. If 
7 is the mean electron Lorentz factor, this occurs at typical energies ~ meC^/7. Assuming 
mildly relativistic particles, the previous limit implies that conservative scattering should 
provide good accuracy up to photon energies of some tens of keV. This, together with the 
fact that resonant scattering occurs in regions where B ^ Bcrit-, makes the use of the (much 
simpler) non-relativistic (Thomson) cross section adequate. We anticipate here that, albeit 
supported by physical considerations, this provides only a zeroth level description and a 
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more thorough treatment demands for the full QED cross section, as it is discussed in more 
detail later on (see A further simplification arises because, under the typical conditions 
encountered in the magnetosphere, vacuum polarization dominates over plasma effects. In 
this situation, the two (ordinary and extraordinary) normal modes are linearly polarized. 

Since radiative de-excitation occurs on a very short timescale, one can safely assume that 
the electron is initially in the ground state. For a particle initially at rest, the non-relativistic 
scatter ing cross secti ons at resonance are easily derived from the general expression given 
e.g. by iHeroldl (119791 ) by performing the substitution 
11 

(5) 



whe re F = (Ae^ujli) / (SmpC^) accounts for the finite transition lifetime of the excited state 



(e.g. 



Daugherty fc Ventura 



1978 



Ventura 



19791 ). Since in the present case it is Hub — 



keV, the resonance peak is so narrow and prominent that non-resonant contributions to the 
cross section are negligible. One can therefore take the limit 



lim -, tt; TTT- 

r-o (cu - UBy + FV4 

which results in 

SvrroC 



2it6{uj — ub) 



(6) 



da 

dQ' 
da 

dn' 
da 

dn' 
da 

dn' 



S{uj — ojb) cos"' 6 cos^ 6' 
5{uj — ujb) cos^ 6 



1-1 



1-2 



2-2 



2-1 



37rroC 



(7) 



SvrroC 



5{uj ~ lob) 



Snroc 



5{uj — ujb) cos^ 6' 



where 9 [9') is the photon angle before (after) the scattering and tq is the classical electron 
radius. Here and in the following the index 1 (2) stands for the ordinary (extraordinary) 
mode. 

Upon normalization, the previous expressions give the probability that an incident pho- 
ton with polarization state i and direction 9 is scattered at angle 9' with polarization state 
j. The total cross sections for separated processes are easily computed by integrating the 
previous expressions over all outgoing photon angles 



0-1-1 
(71-2 



/ dn' 


da 






dn' 




da 


j dn' 






dn' 



1-1 



1-2 



vrroc 
Snroc 



S{uj — ujb) cos^ 9 
5{iJ — lub) cos^ 9 



(8) 
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SvrroC 



0"2-2 



/ dn' 


da 




JiTT 


dn' 




da 


f dVt' 






dn' 



2-2 



2 

TITqC 



S{uj - LUb) 

6{uj - lob) ■ 



2-1 2 

The total cross section for scattering of an incident ordinary (extraordinary) photon is 
obtained by summing the first (second) pair of expressions in equation ([8]). Finally, in order 
to determine the photon direction after scattering (i.e. the two angles 6', 0') in the Monte 
Carlo code, the following integrals are required 



1 fn'^,,d. 



ai^j Jo Jo dn' 

I 1-2-K re' 

^^'^ 
o"i_i JO Jo ail 

— r r^f'^ 

ai_2 Jo Jo dn' 



1-1 



1-2 



2^ 

(72-1 Jo Jo dn' 
f f dn'^ 



a2-2 Jo Jo dn' 



2-1 



2-2 



-(l-cos^^O 
^(1 -cos^') 



(9) 



2.4 Photon propagation in the magnetosphere 

The scattering cross sections discussed in §2.31 hold in the electron rest frame (ERF). In 
particular, both the photon (uj) and the cyclotron {ujb) frequency entering expressions ([7])- 
(Q are evaluated in the ERF. In the case of a charge moving with velocity v = (3c and 
Lorentz factor 7 with respect to a frame attached to the star, the total cross sections (eq. 
[8]) take the form 



0-1-1 
0"2-2 

where 



3^1-2 

3(T2-1 '- 



LOb 



2 



S{uj — cud) cos^ 9 
6{uj - ujd) 



(10) 



7(1 -/5/i)' ^^^^ 
6 is the angle between the incident photon direction and the particle velocity as measured 

in the ERF and /i is the cosine of the same angle but measured in the stellar frame. The 

latter two quantities are related by the usual transformation 



cos 9 



(12) 



1 ■ 

Since particles are moving along 5, the magnetic field is unaffected by the Lorentz trans- 
formation, and the value of B as measured in the stellar frame can be used to compute 
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the cyclotron frequency tu^ in the ERF. It is worth stressing that uo in eqs. (11 01) is now the 
photon frequency in the stellar frame. 

The scattering optical depth for a photon which travels a distance dl in the magneto- 
sphere is 



dTij = d£ 



dpne{r)-f'\l - PlJ.)aij{u, f, /3)/e(f, 7/?) 



(13) 



where the factor 1 — appears because of the change of reference between the ERF and 
the stellar frame, is the (velocity integrated) particle density and fe is the (normalized) 
momentum distribution as defined in equations ([3]) and (jlj). 

The indices i and j refer to the initial and fin al photon polarization sta t es and [Pmin, Pmax] 
is the charge velocity spread. As pointed out by [Fernandez fc Thompson! (120071 ). the integral 
in equation ( |T3i) can be readily calculated by exploiting the (5-function in the scattering cross 
section. Denoting by 



1,2 



2^, / ,2 /i±— V(W^)^ + /^^-l (14) 

the two roots of the quadratic equation uj = ud, the 5-function in frequency can be trans- 
formed into a (5-function in velocity 



6{uj - ud) 



5{P - Pk] 



E 

fc=1.2 



Accordingly, the total scattering depth can be expressed as 



dri 
and 

dT2 



dri^i + dri_2 



2tt roC — ^( 



(I 



dT2-2 + dT'^ 



2-1 



27r^roC 



n^ujB 



fc=l,2 



'E ■ 

fc=l,2 



'1 



Ia* - Pk 



fe{r,-fkPk) 



fe{r,^kPk) 



(15) 



(16) 



(17) 



for photons initially in the p olarization state 1 and 2 , respe ctively. These expressions are 
analogous to those derived by [Fernandez fc Thompson! (120071 ). although their equation (33) 
seems to contain an error. The spatial distribution of charged particles depends, in fact, 
on their average velocity, i.e. the speed at which charge carriers flow, and not on the velocity 
of the single particle. For this reason rif. is not evaluated at P = Pk, but at (P) which is, 
in general, a function of position (see also §2.21) . Moreover, their equation (13) contains an 
unexpected factor u;^) in place of ub- The reason for this is obscure since the ratio ujb/B 
turns out to be independent of both the magnetic fleld and photon energy. 

Once the initial photon polarization, energy and direction have been flxed, equation ([T6j) . 
or (fT7|) . is integrated along the photon path until a scattering occurs (see §3.2p . Although 
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general relativistic effects are certainly important, here we restrict ourselves to newtonian 
gravity and assume that photons move along straight lines between two successive scatter- 
ings. Proper inclusion of null geodesies in a Schwarzschild space-time, albeit conceptually 
simple, turned out to be computationally quite costly and we decided to dismiss it. As it 
is apparent from equation f|T5|) . resonant scattering may occur only when the roots (3k are 
real, i.e. only if {ubIoj)"^ + /i^ — 1 ^ 0. Since ojb depends (through B) on position alone, at 
every point in the magnetosphere the previous co ndition discriminates those pa irs of photon 



energy and angle for which scattering is possible (iFernandez &: Thompson] 120071 ). In case the 
particle velocity is always of a given sign (charge carriers all positive or negative), only the 
roots (3k with the same sign are meaningful. If there exist two roots with the right sign (i.e. 
both are positive or negative), the criterion for selecting onto which particle (the one with 
velocity (3i or (32) the photon actually scatters is discussed in §3.21 

2.5 Seed photon distribution 

Primary photons are assumed to be emitted by the cooling surface of the neutron star. 
Although, up to now, no detailed model for surface emission from a magnetar has been 
presented, it seems unlikely that the spatial and energy distribution of the surface-emitted 
photons are the same as in ordinary cooling neutron stars. In particular, being the surface 
heated by returning currents (e.g. TLK), the surface temperature is expected to be inho- 
mogeneous (with the equatorial belt hotter than the polar regions) and it is unclear if a 
standard (i.e. in hydrostatic ar id radiative equilibrium) atmo s phere can be present on th e 



top of a magnetar (see however iGiiver. Ozel fc Lyutikov II2006I : iGiiver. Ozel fc Gogusll2007l ) 



On the wake of this, in order to keep our treatment as general as possible, we do not 
prescribe an a priori surface temperature distribution (see §3.11) . In the present version of 
the code, the initial energy distribution is taken to be planckian for either ordinary and 
extraordinary photons, although other spectral distributions can be easily accommodated. 
Different degrees of polarization of the primary spectrum can be then obtained adding 
together, in different proportions, ordinary and extraordinary blackbody photons. Since the 
non-isotropic opacity of the stellar crust might convey radiation in a preferred direction, we 
introduce a beaming parameter 6^1 such that the specific intensity at the star surface 
takes the form 

exY>{hv/kT) 



nM oc ^J^'^ (18) 
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where /i is the cosine of the angle between the initial photon direction and the magnetic 
field. For 6 = 1 the radiation is emitted isotropically in the outward hemisphere. 

3 THE MONTE CARLO METHOD 

The code is structured into four main blocks, as outlined below. In the first thermal photons 
are emitted from the stellar surface, in the second the program evaluates the optical depth 
of the photon as it propagates through the magnetosphere, while the third is finalized to 
solve the kinematics of the electron-photon scattering. Finally, escaping photons are stored. 
Each block is briefiy described in the following. 

3.1 Photon emission 

Because of the intrinsic asymmetry of the model, the observed spectrum depends on both 
the shape, and the (longitudinal) position of the emitting region on the star surface, and the 
viewing direction. Moreover, as mentioned above, the star surface temperature distribution 
may not be isotropic. To account for these effects, the star surface is divided into Nq x Nq, 
zones by means of an equally spaced cos G and $ mesh, where G and $ are the magnetic 
colatitude and longitude. This choice guarantees that all patches have the same area, so 
that the number of emitted photons depends only on the patch temperature (i.e. patches at 
the same temperature emit the same number of photons). A different temperature may be 
attached to each surface patch in such a way to reproduce (up to the accuracy allowed by 
the finite mesh resolution) any kind of thermal surface map. 

Initially we fix the coordinates of an emitting patch and assign a value for the polarization 
state s of each seed photon, i.e. s = 1 for the ordinary mode or s = 2 for the extraordinary 
mode. All photons are emitted at the patch centre P. Then, a photon is extracted at random 
from the distribution ( |T8l) . We assume that the initial photon angles are such that the 
azimuth (as referred to B in P) is uniformly distributed while fi = cos6p is obtained solving 
the equation fi = {UpY, where Up is an uniform deviate and h is the beaming parameter 
introduced in eq. f|T8l) . The coordinates of the emission point and the initial momentum 
univocally determine the ray along which the photon moves. 

Actually, after experiencing scattering(s), some photons will reach the star surface again. 
Their number is fairly limited, since scattering typically occurs at a distance Rgc of a few 
stellar radii. The star disc, as seen from the last scattering point, subtends a solid angle 
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~ [Rns/Rsc)'^ ^0-1) cind this is also an upper limit to the fraction of photons which are 
scattered back onto the star surface. Numerical simulations show that the actual value is 
quite smaller, < 1%. We assume that all photons impinging on the surface are absorbed 
(regardless of their polarization state). 

3.2 Scattering depth 

In a Monte Carlo scheme the distance £ a photon of polarization state s travels between 
two successive interactions (i.e. emission-scattering or scattering-scattering) is estimated by 
integrating the scattering depth dTg given by equations (fTB]l and (fT7|l until 



where U is an uniform deviate. Direct numerical evaluation of the integral (fT9l) proved, 
however, quite time consuming, and we found more efficient and faster to perform a step- 
wise integration the differential equations (|T6l) and ( ITTll using a fourth order Runge-Kutta 
method. Integration is terminated as soon as the value of the optical depth exceeds —InU 
and a linear interpolation between the last two steps is used to determine with better accu- 
racy the value of i where = — In U. 

At each integration step we check if the photon still lies in the region of the {uh/uj, fi) 
plane where resonant scattering is allowed, i.e. if uj^/uj'^ + /i^ — 1 ^0 (see §2.4p . When 
the previous inequality is found to hold no more, we further check if the photon trajectory 
is bound to bring it back into the scattering permitted region or not. This is achieved 
by computing numerically the tangent to the photon path [in the (cu^/cj, fi) plane] where 
dj^/cu^ + — 1 ~ and checking if it lies in between the two limiting values {fi±l)/{ujB/uj)- 
If not, the photon is taken to freely escape to infinity [see also figure 1 of Fernandez & 
Thompson (2007)]. The values of the energy and direction of the photon are then stored, 
the program returns to step 1, and a new seed photon is emitted. 

3.3 The scattering process 

Assuming that equality r<j = — In f/ is verified at some distance £ from the point of the 
previous photon interaction, the kinematics of the scattering must be solved in order to 
obtain the new direction and energy of the photon. This requires the knowledge of the 
velocity j3k of the resonant electron and the new photon polarization state. This is obtained 
by generating two new random numbers, Ui and U2, and comparing them with the ratios 




(19) 
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of the corresponding cross sections. For a photon initially in the ordinary polarization state 
(s = 1), mode switching upon scattering occurs if Ui > (Ti_i/((Ti_i + (T1-2) = 1/4, while for 
an initially extraordinary photon {s = 2) this happens if Ui > cr2-2/(o"2-2 + o"2-i) = 3/4. 
Similarly, the decision about onto which of the two resonant electrons (assuming that both 
values of Pk are acceptable) scattering actually occurs is reached by comparing U2 with the 
ratio Ss{Pi)/[Ss{Pi) + Ss{P2)], where Ss{Pk) stands for each addendum in the sum at left 
hand side of equation f|T6|) (fT7]). If U2 < Ss{l3i)/[Ss{l3i) + Ss{l32)], the scattering electron 
velocity is jSi, otherwise it is l32- At this stage, all parameters entering the differential cross 
section of the process are known. 

Upon scattering with a moving charge, the momentum and energy of the photon are 
modified. Since the cross section ([9]) are defined in the ERF, the evaluation of the scattering 
angles 6' and 0' requires a Lorentz transformation from the stellar frame to the frame 
comoving with the resonant electron jSk- For linearly polarized incoming light the distribution 
of the azimuthal angle is isotropic, so that 0' = 2ttU3, where U3 is an uniform deviate. 
Concerning the scattering angle, we note that in the non relativistic case all quantities ([9]) 
are proportional either to 1 — cos^' or to 1 — cos^^'. Then, after drawing a new uniform 
deviate Ki, the scattering angle is given by cos^^' = 2U4 — 1 or cos^ 9' = 211^ — 1, depending 
on the case. 

The corresponding angles in the stellar frame and, hence the new photon direction, are 
obtained by means of Lorentz transformations. In this frame the photon frequency is given 

by 

u' = -ilu:{l-Pk^^){l + PkCose'). (20) 

Finally, once energy and momentum of the scattered photon are known the computation 
proceeds starting again from point 3.2, integrating equations (|T6l) or f|T7j) along the new 
photon path. 

3.4 Photons storage 

Escaping photons are collected on the "sky at infinity", i.e. on a spherical surface located 
sufficiently farther out to see the star (and its magnetosphere) as point-like. We introduce 
an angular grid (O^, $s) which divides the "sky at infinity" in a fixed number of patches, 
similarly to what has been done for the stellar surface. When the escape condition (see 
§2.4p is met, the two angles B and $ which characterize the ray relative to the star centre 
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are computed from the photon momentum and the sky patch hit determined. Counts are 
stored in a three-dimensional array, the first two indices of which label the sky patch while the 
third the photon energy. This allows to analyze the resulting spectra in different directions of 
observation when a large number of events are processed. Each run involves Nt^t photons, and 
is performed changing the initial polarization states s and the co-ordinates of the emitting 
patch. The resulting spectrum is obtained by superposition of the various emitting patches. 

4 RESULTS 

Our Monte Carlo code, written in FORTRAN90, proved to be efficient and relatively fast. 
Despite the complexity of the whole procedure, we can process about 7000 photons/s on 
a dual-core Xeon 2.8 GHz machine. The CPU time for a typical production run (several 
million photons) is 10-20 min. We stress that the result of each run is a 3D array which 
gives the number of counts at different positions on the sky and at different energies (see 
§3.4p . Further manipulations (e.g. to account for viewing angles, or to derive the pulse shape, 
see below) are performed at the post-processing level by means of IDL scripts, at negligible 
computational cost. In the following subsections we discuss the general properties of our 
spectral models. 

4.1 Spectra 

In order to explore the role of the different parameters we computed a set of spectra, by 
evolving Npatch = 150, 000 photons for A^e x = 8 x 4 = 32 surface patches (i.e. each model 
has Nfot = 4, 800, 000 photons). We assume that the star surface is at constant temperature, 
and that the seed radiation is isotropic (6 = 1, see § 12.51) and completely polarized, either in 
the ordinary or extraordinary mode. Furthermore, we treat the case of an aligned rotator, 
i.e. the spin and magnetic axes coincide. Photons are collected onto a A^^e^ x A^$s = 10 x 10 
angular grid on the sky, and in Ne = 50 energy bins in the range 0.1-100 keV. The magnetic 
field has been fixed at B = 10^^ G and the surface temperature at kT = 0.5 keV. The mean 
and the maximum number of scatterings per photon are in the ranges ~ 0.5-2 and ~ 10-20, 
respectively, depending on the parameter values and on the location of the emitting patch 
on the star surface. 

In Fig. [1] we show the spectra, averaged over as seen by observers whose line-of-sight 
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(LOS) is at different angles 0^ with the star spin axii^. The most sahent characteristic is the 
absence of symmetry between the north and the south hemispheres: as increases, spectra 
become more and more comptonized. This reflects our choice for the electron velocity distri- 
bution, which accounts for the charges bulk velocity, and currents flow from the north to the 
south pole along the field lines (of course the opposite choice for the current direction would 
simply result in 0^ 180° — O^). We found that the spectral shape is almost insensitive to 
the seed photons polarization state (see Fig. [1]). This means that observations of the phase 
averaged spectrum are not expected to provide useful insights into the polarization degree 
of the surface emission (but see §4.2p . 

Figs. El [31 m illustrate the effects on the spectral shape of varying Pbuik, kTe and A^, 
respectively (here and in the following we put A(f) = A^jy-s to simplify the notation). 
Spectra have been averaged over $s, and plotted for two values of B^, one for each hemisphere 
(left and right panels). As it can be seen, an increase in each of these parameters (either 
Pbuik, kTe or A(f)) always corresponds to an increase in the comptonization degree of the 
spectrum. The effect is particularly notable in the case of Pbuik- If Pbuik ^ 0.5 an observer 
located in the southern hemisphere (i.e. with currents flowing towards him) sees a spectrum 
which is no more peaked at ~ kT, but peaks instead at about the thermal energy of the 
scattering particles. This is because comptonization starts to saturate and photons fills the 
Wien peak of the Bose-Einstein distribution. For intermediate values of the parameters, 
spectra can be double humped, with a downturn between the two humps (a clear example 
of this behaviour is illustrated in Fig. [5l). We note that some of the model spectra presented 
by Fernandez & Thompson (2007) also exhibit a downward break in the tens of keV range. 
In particular, when assuming a (non-thermal) top-hat or a broadband velocity distribution 
for the magnetospheric charges, they found that multiple peaks can appear in the spectrum. 
The difference is that our model predicts at most two peaks, and that the energy of the 
second on e gives a direct informat i on on the energy of the magn etospheric particles. As 



noticed by [Fernandez fc Thompson! (120071 ) and lEsposito et al.l (120071 ). double peaked spectra 



may play a role in the interpretation of the broadband X-ray spectrum of SGR 1900+14 
and SGR 1806-20. In particular, the detection of a spectral break at about a few tens 
of keV may have remarkable physical implications and provide important diagnostics for 



^ Note that the total number of collected photons is usually lower than Ntot (4,800,000 in the present case) since a (small) 
fraction of photons reach infinity with an energy outside our range of collection (i.e. 0.1-100 keV). 
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the physical parameters of the model. A spectral break at ~ 15 keV, as the one possibly 
detected in the case of SGR 1 806-20, would transla te then in a temperature of ~ 5 keV for 
the magnetospheric electrons f Esposito et al.l 2007 ). 

The efficiency of the resonant scattering also increases by increasing kTe (Fig. [3]), al- 
though this effect is less pronounced than that observed while increasing the current bulk 
velocity. This is expected, because a change in Te corresponds to a change in the average 
thermal velocity for the magnetospheric particles, and not to a boost that equally affects each 
single particle. Similarly goes for A0, which effect is less pronounced than that of the bulk 
velocity (see Fig. H]). Again, we find that no significant spectral change occurs exchanging 
the polarization of the seed photons from ordinary to extraordinary. 

Although it would be inappropriate to define the RCS spectra as a "blackbody plus 
power-law" (the double-humped spectra shown in Fig. [5] are definitely far away from such a 
definition), in many cases the general shape of the continuum is that of a thermal bump and 
a high-energy tail. In this sense model spectra are reminiscent of the empirical blackbody 
plus power-law model often used to fit (rather successfully) the magnetars soft X-ray emis- 
sion. Since, when present, the high-energy tail is indeed power- law- like, it is of interest to 
investigate how the spectral index F (as derived by fitting the high-energy tail with a power- 
law) changes with the parameters. In particular, a hardening of the spectrum is expected for 



increasing twist angle (TLK) and this was invoked as a 
corre lated fiux-hardening variations in some sources (e.g. 



jossible mechanism to explain the 



Mereghetti et al 



2005b: 



-lea et al. 



20051 ). This is confirmed by our calculations (see also [Fernandez fc Thompson 



20071), as 



shown in Fig. El The photon index monotonically decreases with A0, going, in the present 
case, from ~ 3 to ~ 2.4 by changing the twist angle by ~ 1 rad. The behaviour is quite similar 
at both the field strengths we considered, although spectra for B = 10^^ G are fractionally 
harder. The model shown here has kT^, = 30 keV, l3f,uik = 0.3, a uniform surface temperature 
kT = 0.5 keV and spectra have been obtained summing over all the sky patches. 

To illustrate the effects of a non-homogeneous surface temperature distribution, we dis- 
cuss the case in which photons are emitted by a single patch. The subdivision of the star 
surface and of the sky is the same as that adopted before, and also the energy range and 
bin width. In the present run the seed radiation is taken to be isotropic {b = 1) and unpo- 
larized, i.e. an equal number of ordinary or extraordinary photons are emitted, and, again, 
the spin and magnetic axes coincide. We selected an emitting patch located just above the 
equator (centred at = 77.5°, $ = 45°) with a surface temperature of kT = 0.5 keV. The 
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magnetospheric parameters are Pbuik = 0.3, kT^ = 10 keV and A0 = 1.3. Figure [7] shows 
the emerging spectrum, as viewed by an observer whose LOS makes an angle = 90° with 
the spin axis (i.e. the star is seen equator-on) for different values of the observing longitude, 
$s = 20°, 140°, 220°. These three values correspond to having the emitting patch in full 
view (seen nearly face on), partially in view and screened by the star. The effects of the 
different viewing angle on the spectrum are dramatic. When the emitting patch is in full 
view both the primary, soft photons and those which undergo repeated resonant scattering 
reach the observer and the spectrum is qualitatively similar to those presented earlier on, 
with a thermal component and an extended power-law-like tail. On the other hand, if the 
emitting region is not directly visible, no contribution from the primary blackbody photons 
is present. The spectrum, which is made up only by those photons which after scattering 
propagate "backwards", is depressed and has a much more distinct non-thermal shape. 



4.2 Polarization of the emitted radiation 

Radiation emerging from strongly magnetized neutron stars is expected to be highly polar- 
ized, due to the strong dependence of radiation transport on the photon propagation mode. 
Polarization studies have already started at low energies (IR), and future X- and 7-ray 
polarimetry with high sensitivity instruments, such as the planned photoelectric polarime- 
ter to be flown on the ESA mission XEUS, are expected to extend them over a broader 
spectral band. The development of detailed theoretical predictions is therefore fundamental: 
polarimetry will bring into view a new and unique dimension of the problem, through the 
knowledge of polarization degree and swing angle. 

In our scenario, the degree of polarization in the soft X-ray radiation emitted by mag- 
netars results from a combination of several effects. Seed thermal photons, origi nating from 



the crust or atmosphere of the star, do posses an intrinsic polarization (e.g. 



Zane et al 



2000 



Van Adelsberg fc Laill2006l ). The fraction of polarization, which is determined by the 
competition between plasma and vacuum properties, depends on the energy band, and on 
the details of the density and temperature gradient in the emitting region. Seed photons 
then propagate in the magnetosphere, where multiple resonant scatterings further influence 
the polarization degree. By using our Monte Carlo simulation, we are in the position to 
investigate the latter effect, i.e. to estimate the degree of polarization which is expected to 



20 L. Nobili, R. Turolla and S. Zane 

arise because of magnetospheric effects only and to investigate its dependence on the model 
parameters. 

In Figs. [H] and [9] we show, as a function of various parameters, the degree of polarization 
of the emerging radiation, defined as \Nextr — Nordl/iN^xtr + ^ord) where N^xtr and Nord 
are, respectively, the number of ordinary and extraordinary photons collected at infinity. 
The polarization degree has been averaged over frequency, over the whole emitting surface 
and over the sky at infinity. As it can be seen, the efficiency at which completely polar- 
ized surface radiation is de-polarized increases by increasing the strength of magnetospheric 
upscattering, i.e. by increasing one of the three parameters Pbuik, kTe or A0. This effect is 
stronger for ordinary seed photons, for which the probability of undergoing mode switching 
in the scattering process is higher (see, e.g., eqs. [10]) and for photons emitted close to the 
south pole (see Fig. [10], the latter result reflects our choice for the direction of the current 
flow, as discussed earlier). On the other hand, would the surface radiation be completely 
unpolarized, we can see that, while passing through the magnetosphere, it can acquire only a 
relatively small degree of linear polarization: typically 10-20%, up to 30% for very extreme 
values of the current bulk velocity. This means that, would future observations of X-ray 
polarization result in measurements larger than 10-30%, the excess has to be attributed to 
an intrinsic property of the surface radiation. 

We have also explored how the polarization degree depends on the photon energy and 
a representative case is shown in Fig. [HI The two panels refer to a run with the same set 
of model parameters {Ptuik = 0.3, kTf, = lOkeV, kT = 0.5 keV, A0 = 1.3) but performed 
assuming that seed photons are completely polarized either in the extraordinary (left panel) 
or ordinary mode (right panel). The polarization degree has been computed as above, but 
now different viewing directions are retained (i.e. only sum over $s has been performed). 
Emission is again from the entire surface (at constant T) and the star is an aligned rotator. As 
expected, for 100% polarized seed photons the polarization degree decreases with increasing 
energy, since harder photons undergo more scatterings. Low energy photons tend to keep 
their original polarization state, although there is a dependence on the viewing angle. Not 
surprisingly, even at low energies, the polarization degree is higher when the LOS is close to 
the north pole (dash-triple dotted lines in Fig. [TTl) and drops for increasing viewing angle. It 
is interesting to note that the largest de-polarization (at low energies) does not occur close 
to 180° but when viewing the star southern hemisphere at an intermediate angle because of 
the low particle density near the poles. 



X-ray spectra from magnetar candidates 21 

4.3 Viewing angle effects 

Spectra presented in §4. II have been computed accounting for different viewing angles only in 
the case in which the star is an aligned rotator, i.e assuming that the spin and magnetic axes 
coincide. Under this hypothesis, the viewing geometry is described by a single angle which is 
just the colatitude of the centre of the sky patch where photons are collected. Since the 
magnetic field and the current distribution are axially symmetric, the contributions from all 
the sky patches located at the same value of (and different $s) may be summed together 
if surface emission is itself axisymmetric, as in the uniform temperature case discussed at 
the beginning of this section. 

In order to treat the more general case in which the spin and magnetic axes are not 
aligned, we introduce two angles, x ^^id ^, which give, respectively, the inclination of the 
LOS and of the dipole axis with respect to the star spin axis. This also allows us to take into 
account for the star rotation and hence derive pulse shapes and phase-resolved spectroscopy. 
Because of the lack of north-south symmetry, it is ^ x ^ while ^ spans the interval 
[0, tt/2]. By introducing the rotational phase a (0 ^ a ^ 2tt), the co-ordinates of the point 
which represents the intersection of the LOS with the sky for each value of a are 



cosBs = cos X cos ^ + sin X sin ^ cos a (21) 
— cos X — cos Qs cos ^ 

cos $s = = . 

sin Bs sin ^ 

At constant x ^i-nd ^, eqs. (1211) trace a circle on the sphere which represents the sky. As 
a result of each Monte Carlo run, the spectrum in counts has been recorded for each pair of 
values Qg^i, which correspond to the centres of the sky patches, N{Qs,i, ^sj, Ek). In order 
to compute the spectrum at a discrete set of phases ai, we perform a double interpolation of 
this array over the angular variables, to obtain the number of counts in correspondence to 
the pair of angles Qs{ai), ^s{c(i) given by eq. ([21]), i.e. Nph{ai, Ek). Finally, integration of Nph 
over E or a gives the lightcurve in a given energy band, or the phase-averaged spectrum, 
respectively. An illustration of the effects of a different viewing geometry is shown in in 
Fig. [121 where spectra correspond to increasing values of x- 

A systematic investigation of the properties of the pulse shape while varying the model 
parameters is beyond the purpose of the present paper, and it will be presented elsewhere 
(Albano et al. in prep.). Here we just show in Fig. [13] two examples, both relative to a star 
seen equator-on {x = 90°), but for two different inclinations of the magnetic axis (^ = 10° 
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and ^ = 50°). In the first case the pulse profiles in the soft (0.5-2 keV) and hard (2-6 keV) 
band are shifted in phase by ~ 180°. By increasing ^ the pulsed fraction and the pulse shape 
sensibly change with the energy band. The pulsed fraction increases with the energy and, 
at the same time, the double peaked structure present in the low energy band disappears at 
higher energy where the lightcurve is sinusoidal. 



5 XSPEC IMPLEMENTATION AND APPLICATIONS 

One of the goals of the present investigation is to apply the resonant compton scattering 
model discussed in the previous sections to magnetar spectral fitting, by implementing it 
into the standard package for X-ray spectral data analysis XSPEC. Clearly, our Monte Carlo 
spectra can be loaded in XSPEC only in tabular form, using the at able option. This implies 
that a model archive has to be generated beforehand, for a reasonably wide range of the 
model parameters. Although a production run takes (under typical conditions) about 20 
m, building a large model archive necessary demands for a compromise between generality, 
accuracy and feasibility. As we discussed already (see §4.11) . the model has four parameters: 
(3buiki kTe, A(j) and kT, assuming that the surface is at constant temperature. If a model 
is computed for, say, ten values of each parameter, this would result in a total of 10^ runs 
requiring about 2 x 10^ m ~ 140 d of CPU time. Even splitting the computation over a 
few machines, the time needed (~ month) is barely acceptable. Moreover, we are aware that 
the adopted description of the charge velocity distribution, which involves two out of four 
model parameters, is far from being consistent. For these reasons, we decided to simplify our 
treatment by imposing that the electron bulk kinetic and thermal energies are related. The 
mean thermal energy for a ID relativistic Maxwellian distribution can not be expressed in 
closed form. However, to an excellent accuracy, it is 

(7-i>-^i7|ky (22) 

We then derive the value of the electron temperature by assuming equipartition between 
thermal and bulk kinetic energy, i.e. by solving for Og = /cTe/mgC^ the equation 

Ibuik - 1 = 21/(1+6,) • (23) 

In order to avoid that for the higher values of jSbuik we consider (see below) the assumption of 
conservative scattering in the ERF is invalidated, the solution of eq. fl23l) is actually halved. 
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The grid of models has been generated for ^ A0 ^ 2 (step 0.1), 0.1 ^ Pbuik ^ 0.9 
(step 0.1) and eight values of kT (in keV): 0.1, 0.13, 0.16, 0.2, 0.25, 0.40, 0.63, 1, under the 
assumptions that the surface has a constant temperature, emits isotropically and the surface 
radiation is unpolarized. The number of divisions on the star surface and on the sky, the 
energy range and bins are taken as in § I4.H but now we evolve Npatch = 225, 000 photons per 
patch, therefore each model corresponds to Ntot = 7, 200, 000 photons. Again, the magnetic 
field is fixed at B = 10^^ G (further archives corresponding to different values of B can be 
easily generated). The computed spectra have then been averaged over the whole sky at 
infinity, smoothed and re-interpolated (using a logarithmic interpolation) over a grid of 300 
equally spaced energies in the range 0.1-15 keV and on a logarithmic grid of 100 equally 
spaced temperatures in the range —1 ^ log/cT ^ 0. The latter step is necessary because 
interpolation on the logarithm of the spectrum with respect to parameters is not possible 
within XSPEC for tabular models. After some experimenting, we found that in order to have 
enough accuracy when interpolating the spectrum a fine grid in kT is necessary. The final 
XSPEC atable spectral model (22 MB in size, named ntznoang.mod) has been created by 
using the routine wf tbmd, available on-line|^ 

The ntznoang model has four free parameters {(3buik, logkT plus a normalization 
constant), which can be simultaneously varied during the spectral fitting following the stan- 
dard minimization technique. It is important to note that this model has the same number 
of free parameters than t he canonical blackb o dy plus power-law empirical model or the RCS 



model recently discussed in lRea et al.l (120081 ). and hence has the same statistical significance. 

Following essentially the same procedure outlined above and making use of the same 
archive, we have also built a XSPEC model in which the dependence on the two geometrical 
angles, x and C,, is explicitly accounted for, as discussed in §4.31 Phase-averaged spectra have 
been computed on a 7 x 7 equally-spaced grid of x and ^ values. The two angles are in the in 
the ranges ^ x ^ 180° and ^ ^ 90°, respectively. At variance with the angle-averaged 
case considered previously, the grids in the other parameters (except kT) are coarser: ^ 
A0 ^1.8 (10 values, step 0.2), and 0.1 ^ Pbuik ^ 0.9 (5 values, step 0.2). Maintaining the 
same parameter grids used to build the ntznoang model would, in fact, result in too a large 
file to be read into XSPEC. The final atable spectral model, ntzang.mod, is ~ 300 MB 



^ see |http://heas arc. gsfc. nasa.gov/docs/heasarc/ofwg/docs/gener al/1 



modelfiles_memo/modelfiles_memo.html. 
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in size, and has six free parameters {(3buik, ^4', log/cT, x, ^ plus a normalization constant). 
Despite the larger number of free parameters, the "angular" model can be used to infer 
information about the viewing geometry, eventually combining information that can be 
obtained by fitting simultaneously phase-resolved spectra, or independently from the study 
of the pulse profile. 



A systematic application of both models to magnetars spectra is in progress, and will 
be reported elsewhere (Israel et al. and Rea et al. in preparation). Here we present only an 
example which is illustrative of how the two atable spectral models behave when applied 
to X-ray data. Fig. [H (left panel) shows the fit of the 0.1-10 keV XMM-Newton EPIC-pn 
spectrum of the transient AXP CXOU J1647-4552 taken on February 17 2007, i.e. about 



five months after a burst and a glitch were detected from th is source (IKrimm et al 



2006; 



Israel et al. 



2006. see also 



Muno et al. 



2007 



Israel et al. 



20071 ). All details about the obser- 



vation will be reported in Israel et al. (in preparation). The spectrum has been modeled with 
the angle-integrated ntznoang model, modified by interstellar absorption (phabs model in 
XSPEC). Data and best fitting model are shown in Fig. [13] and the best fit parameters are 
listed i n Tabled As exp ected, and as it has been also found in other applications of the RCS 



model ( Rea et al. 



20081 ) , the inferred value of the column density is smaller than that implied 



by a blackbody plus power-law fit. This is because the empirical blackbody plus power-law 
modelling is known to overestimate the soft X-ray emission and, in turn, the value of the 
interstellar absorption. 



Since the fit is already very good {xl = 0.81), there is no statistical need to introduce 
two further parameters. However, we also tried to fit the same observation with an absorbed 
ntzang model, with the only goal to check and test the correctness of its XSPEC imple- 
mentation; results are shown in Fig [H] (right panel) and reported in Table [H As expected, 
the values of the angles are unconstrained, and the remaining parameters are in agreement 
with those found with the first model. Again, is not our main scope to provide the physical 
values of the angles here: instead we stress that this figure is presented purely as an illustra- 
tion. Nevertheless, the successful spectral fit with the ntznoang model clearly demonstrates 
that the model can catch the main features of the magnetar emission and reproduce them 
quantitatively. 
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6 DISCUSSION AND CONCLUSIONS 



In this paper we have investigated how the thermal spectrum emitted by the star surface gets 
distorted by repeated resonant scatterings onto mildly relativistic magnetospheric electrons 
using a Monte Carlo technique. The goal of this study has been twofold. Our first motivation 
has been to create a model archive which could be implemented as a tabulated model 
in XSPEC and directly applied to fit the spectra of magnetar candidates. The model is 
available in two versions, with or without the explicit dependence on the two angles which 
give the inclination of the line-of-sight and the magnetic axis wrt the star spin. A systematic 
application to different sources is under way and here a (preliminary) fit to the XMM-Newton 
spectrum of CXOU J1647-4552has been presented, mainly for illustrative purposes. 

In building our Monte Carlo c ode we have followed an approach similar to that discussed 



m 



Fernandez fc Thompson! (120071 ). However, the two codes differ in many respects. A major 



difference is in the adopted description of the velocity distribution of the scattering par- 
ticles. We have explicitly accounted for the collective (bulk) electron motion associated to 
the charge flow in the magnetosphere, superimposed to which we assume a ID relativistic 
Maxwellian distribution which simulates the particle velocity spread. We also allow for a 
completely general description of the star surface thermal map and this makes it possible 
to assess the effects of a (spatially) localized emission (e.g. by a hot spot). Moreover, in our 
treatment seed photons are not taken to move only in the radial direction but are drawn 
from a prescribed angular distribution which can account for magnetic beaming effects. 

A s the present app l icatio n to CXOU J1647-4552 shows (§5] see also iLyutikov &: Gavriil 



2006 



Rea et al. 



2007a. 



20081 ). spectral models based on resonant cyclotron up-scattering of 



thermal photons in the magnetosphere of magnetars prove quite successful in interpreting 
quantitatively the soft (~ 1-10 keV) emission from AXPs and SGRs. Albeit the numerical 
computation presented here includes several important details about the microphysics and 
the magnetospheric properties and geometry, it relies on some simplifying assumptions which 
reflect our poor knowledge on some key issues of magnetar physics. 

A prominent one is the nature of the plasma which fills the magnetosphere. Most investi- 
gations on RCS, including our, restricted to unidirectional flows, i.e. assumed that scattering 
occur onto electrons ( a simple bi-directional flow was co nsidered by lFernandez fc Thompson 



20071 ). As discussed by lBeloborodov fc Thompson! (12007! ). in a twisted magnetosphere charges, 
accelerated by the self-induction electric field, may produce e^. Pairs definitely contribute 
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to the scattering depth [f]. The final spectral shape depends on which species populate the 
corona and on their spatial and velocity distribution. Our choice of modelling the e~ current 
in terms of a bulk motion plus a velocity spread seems to be at leas t in qualitative agreement 
with the analysis presented by iBeloborodov fc Thompson! (120071 ) . We point out, however, 
that the assumption of a ID thermal distribution for the particle velocity in the local rest 
frame is somehow arbitrary and no attempt has been made here to assess the effects of 
other possible (local) di s tribut ions. This has been done, in a few representative cases, by 



Fernandez fc Thompson! (!2007! ). who did not include, however, the charge bulk motion. By 



comparing our results with their, one may conclude that, while the general effects induced 
by magnetospheric RCS on primary thermal photons (i.e. the formation of a "thermal-plus- 
power-law" spectrum) are not much sensitive to the assumed particle velocity distribution, 
the details of the spectral shape do. 

A further caveat concerns the star temperature distribution and the primary spectrum. 
Our model archive has been generated assuming that the star radiates a blackbody from a 
uniformly heated surface. At present it is unclear if magnetars do possess an atmosphere. 
A possibility is that highly energetic electrons hitting the surfac e knock out protons which 
then sublimate giving rise to a "current induced" atmosphere ( IBeloborodov fc Thompson 



200?! ). Departures from a blackbody primary spectrum due to reprocess ing in a strongly 



magnetized atinosphere are, how ever, not expected to be dramatic (see e.g. 



Zane et al 



2001 



Ho fc Lai 



2001 



Lai fc Ho 



2003! ) . On the other hand, the issue of the surface thermal map 



appears more serious since even passively cooling iso 



non-uniform surface temperature (see e.g. !Page 



1995 



ated neutron stars are known to have a 



Zane fc Turolla 



2006! ). In the case of a 



magnetar, returning currents impacting on the star surface produce localized heating (TLK). 
Moreover, starquakes, possibly triggered by the strain accumulated during t he growth of the 



twist and connected to the glitching activity discovered in AXPs (see e.g. !DairOsso et al. 



2003! ). can further contribute to the injection of heat into limited portions of the crust. 



Transient AXPs might be powered in a similar way by the sudden release of energy into 



As discussed by 



Medin fc Lail 1 20071 ) . for an iron crust and magnetic fields as fiigh as ~ lO^^G, vacuum gaps may be formed 



above tfic polar regions of SGRs/AXPs, with subsequent pair creation. The pair-dominated region, however, is very thin and 
located just above the star surface. This implies that scattering is resonant for photon energies in the tens of MeV range. 
Since thermal emission from the star surface does not supply such high-energy photons, pair cascades produced by the gap 
breakdown are not going to affect our results. 
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a localized area of the star surface, as observations of the TAXP XTE J1810-197 seem to 



indicate (IGotthelf fc HalpernI 120071 ). 

Although the twisted dipole model used here has the advantage of simplicity while catch- 
ing the essential physical features, most probably it gives only an idealized representation of 
the magnetic field outside a magnetar. The twist may be confined at high magnetic latitudes 
(TLK), or, if global, it might involve magnetic configurations more complex than a dipole. 



Possible evidence for a twist which inv olves in the first p 



magnetic poles have been discussed by 



Woods et al. 



ace the field lines closer to the 



(120071 ) in connection with the period 



derivative evolution and its correlation with spectral hardness in SGR 1806-20 before and 
after the giant fiare of December 2 7 2004^ 



Both 



Lyutikov fc Gavriill (120061 ) and [Fernandez fc ThompsonI (120071 ) assumed that scat- 



tering is conservative in the electron rest frame. As discussed in §2.31 this choice is quite 
adequate if spectral modeling is restricted to the soft X-ray range and has been retained in 
the present work. However, the X-ray spectra of magnetar candidates are nowadays known 
to exhibit also a high energy (~ 20-200 keV) component, which is completely non-thermal 
and is responsible for about half of the bolometric fiux. Although different scenarios for the 
origin of the high energy emission from magnetars have been put forward, not necessary 
involving RCS (see ^II), an intriguing possibilit y is that also the hard tail arises because 
of resonant upscattering in the magnetosphere (IBaring fc Harding! 120071 ) . Given the much 
higher photon energies (in the 100 keV range) this necessary requires the presence of highly 
relativistic electrons (pairs), and, consequently, any attempt to model RCS under those 
conditions demands for a fully relativistic, QED treatment of the scattering cross sections. 
Although we presented here spectra extending up to 100 keV, they must be considered 
as trustworthy only until huo <^ meC^/7, i.e. up to a few tens of keV. Above these en- 
ergies electron recoil starts to become important and the spectrum is expected to break. 
The precise localization of the break would come only from a consistent treatment, and is 
particul arly important to explain the COMP TEL upper limits observed in some magnetar 
sources ( iKuiper et al.ll2006l : iRea et al.ll2007bl ). Moreover, if hard tails are due to a secondary 
po pulation of ultra-relativi stic electrons confined close to the stellar surface (as proposed 
by iBaring fc Harding! 120071 ) . resonant scattering would occur at much higher values of the 
magnetic field, B > Bqed, which makes the need of a completely QED treatment of the 
cross section even more necessary. Same holds for computations aimed at assessing the role 
of ions in shaping the spectra. As previously discussed see §2.2[ positively charged ions are 
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expected to populate the twisted magnetosphere, but whether these particles can effectively 
shape the X-ray spectra is mainly related to the role of those ions located close to the 
star surface. The inclusion of this effect, however, requires the knowledge of the full QED 
resonant cross section for protons/ions which at present has not been investigated in detail. 

Future work needs to address this issue, among others. Clearly, in order to to include the 
relativistic treatment of the scattering process in the electron rest frame, having a tested, 
reliable Monte Carlo code which can be easily generalized is of fundamental importance 
and this has been our second motivation in undertaking this study. In order to extend 
our computation of resonant electron cyclotron scattering to the relativistic regime, we are 
already completing a detailed investigation of the QED resonant cross section. This will be 
then implemented in our Monte Carlo code and results will be presented in forthcoming 
papers (Nobili, Turolla & Zane in preparation). 
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Table 1. Best fit values of the spectral parameters. 



Parameters 


ntznoang 


ntzang 


Nh 
kT 

Pbulk 
A</> 
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n 625+"-""^ 

U.DZO_„ ggg 

60 +" "^ 
'^■"'^ -0.02 

40+" °=^ 

U-*U_g 32 


'"-0.01 
n cQ +0.07 

0.63 to.oi 
65 

'-'•"0 -0.07 

47 +0 03 

U-^' -0.06 


X 




2.1 ± 1.8 


? 

Norm 


U-08i_g 




Flux 


6 


6 


Xi (dof) 


0.81 (145) 


0.83 (143) 



Errors in the parameters arc at Icr confidence level, Nh is in units of 10^^ cm ^, kT is in keV, Xi^ degrees and the 

observed fiux (1-10 keV) is in units of 10~^^ ergscm"^. 




Figure 1. Left; Computed spectra for B = lO^"* G, kT = 0.5 kcV, kTe = 30 keV, iSttdk = 0.3, A4> = 1 and different values of 
the sky colatitude 6.,: 27° (long dashed), 64° (dasli-triplc dotted), 90° (dash-dotted), 116° (short dashed) and 153° (dotted). 
The solid line rcijrcsents the seed blackbody and counts have been summed over <l>s. Here seed photons arc assumed to be 
completely polarized in the ordinary mode. Right: Same, but for seed photons completely polarized in the extraordinary mode. 



This paper has been typeset from a T^/ E^T^ file prepared by the author. 



32 L. Nohili, R. Turolla and S. Zane 




1 2 

log K (keV) 



kT = 0.50 keV 
kT = 30.00 keV 



1 2 

log E (keV) 




5 4 



kT = 50 keV 
kT = 30.00 keV 



log E (keV) 



Figure 2. Top: Computed spectra for B = 10^"- G, kT = 0.5 keV, kTe = 30 kcV, Ac/) = 1 and different values of fft^ik: 0.3 
(dotted), 0.5 (short dashed), 0.7 (dash-dotted) and 0.9 (dash-triple dotted). The solid line represents the seed blackbody and 
counts have been summed over 3>s. The two panels correspond to two different values of the magnetic colatitude: 0s = 64° 
(left) and @s = 116° (right). Seed photons are assumed to be 100% polarized in the ordinary mode. Bottom: Same, but for 
seed photons 100% polarized in the extraordinary mode. 
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Figure 3. Top: Computed spectra for B = 10^* G, kT = 0.5 kcV. 0t„ik = 0.3, A0 = 1 and different values of kTe: 5 keV 
(dotted), 15 keV (short dashed), 30 keV (dash-dotted), 60 kcV (dasli-triplc dotted) and 120 keV (long dashed). The solid line 
represents the seed blackbody and counts have been summed over <i>s. The two panels correspond to two different values of the 
magnetic colatitude: ©s = 64° (left) and ©s = 116° (right). Seed photons are assumed to be 100% polarized in the ordinary 
mode. Bottom: Same, but for seed photons 100% polarized in the extraordinary mode. 
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Figure 4. Top: Computed spectra for B = 10^"' G, kT = 0.5 keV, kTc = 30 keV, I3i,^ik = 0.3 and different values of A(j): 0.3 
(dotted), 0.5 (sliort dashed), 0.7 (dasli-dotted), 0.9 (dasii-triple dotted), 1.1 (long dashed) and Acp = 1.2 (solid line, top). The 
solid line at the bottom represents the seed blackbody and counts have been summed over $s. The two panels correspond to 
two different values of the magnetic colatitude: 0s = 64° (left) and 0s = 116° (right). Seed photons are assumed to be 100% 
polarized in the ordinary mode. Bottom: Same, but for seed photons 100% polarized in the extraordinary mode. 
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Figure 5. Computed spectra for B = lO^'* G, Ac/) = 2; the star is an aligned rotator seen north pole-on. Solid line: kT = 0.1 keV, 
Pbulk = 0-7; dashed line: kT = 0.6 keV, Phulk = 0-6. In both cases fcTe is related to Pbulk through eq. II23I I: seed photons are 
unpolarizcd. 
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Figure 6. Photon index versus A4> for B = lO" G (triangles) and B = 10^^ G (diamonds). See text for details 
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Figure 7. Spectrum from a single emitting patch on the star surface. The LOS is at ©s = 90° and $s = 20° (dotted line), 
140° (dashed line) and 220° (dash-dotted line). The solid line represents the seed bla<;kbody. Because photons are collected in 
a single patch on the sky, the counting statistics is low at the higher energies and the spectrum looks "noisy" . 
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Figure 8. Left: degree of polarization as a function of A<p for B = 10^^ G, kT = 0.5 keV, fibulk = 0.3 and kTe = 30 keV. 
Right: Same, but as a function of kTg for = 1. In both panels different curves correspond to: seed photons 100% polarized 
in the ordinary (solid line), extraordinary mode (dotted line), and unpolarized (dashed line). See text for details. 
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Figure 9. Degree of polarization as a function of /3(,„jj, for kT^ = 30 keV and A<^ = 1. Other parameters and line code as in 

fig. El 




Figure 10. Right: Degree of polarization as a function of the magnetic colatitude © at which seed photons are emitted. Here 
kT = 0.5 keV, /3f,„ifc = 0.3, fceTe = 30 keV and A<f> = 1. Left panel: B = 10^"' G; right panel: B = 10^^ G. In both panels 
different curves correspond to: ordinary (solid line), extraordinary seed photons (dotted line) an unpolarized seed photons 
(dashed line). Photons have been integrated over the entire sky at infinity and over the azimuthal angle at the star surface. 




Figure 11. Degree of polarization as a function of energy for different values of the LOS inclination, @s = 162° (solid), 
126° (dotted), 90° (dashed), 54° (dash-dotted) and 18° (dash-triple dotted). Left: 100% extraordinary polarized seed photons. 
Right: 100% ordinary polarized seed photons. As in Fig. [3 the low statistics is responsible for the noisy appearance of the plot 
at higher energies. The decrease in the polarization degree with energy is clearly visible notwithstanding. 
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Figure 12. Computed spectra for B = lO" G, kT = 0.5 keV, kT^ = 30 keV, /Strife = 0.3, A0 = 1, ^ = 45° and five different 
values of the viewing angle x- 0.01° (dotted), 45° (short dashed), 90° (dash-dotted), 135° (dashed-triple dotted), and 180° 
(long dashed). The solid line represents the seed blackbody. Here seed photons are assumed to be completely polarized in 
extraordinary mode. 
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Figure 13. Left: the lightcurves in the soft (0.5-2 keV, solid line) and hard X-ray band (2-6 keV, dashed line); the model 
parameters are x = 90°, ^ = 10°, Ai/) = 0.7, Pb-ulk = 0-3 and kT = 0.3 keV. Right: same as in the left panel, but for J = 50°. 
See text for details. 
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Figure 14. Left: Fit of the XMM-Newton EPIC-pn spectrum of CXOU J1647-4552 with an absorbed ntznoang model. Top: 
data and best fit model; bottom: residuals. Right: the same observation fitted with an absorbed ntzang model. 



